knitr::opts_chunk$set(echo = TRUE)

rm(list = ls())

#cmdstanr::set_cmdstan_path(path = "C:/Users/kueng/.cmdstan/cmdstan-2.35.0")
#cmdstanr::set_cmdstan_path(path = "C:/Users/pascku/.cmdstan/cmdstan-2.36.0")

library(tidyverse)
library(R.utils)
library(wbCorr)
library(readxl)
library(kableExtra)
library(brms)
library(bayesplot)
library(see)
library(beepr)
library(DHARMa)
library(digest)



source(file.path('Functions', 'ReportModels.R'))
source(file.path('Functions', 'PrettyTables.R'))
source(file.path('Functions', 'ReportMeasures.R'))
source(file.path('Functions', 'PrepareData.R'))

report_function_hash <- digest::digest(summarize_brms)
system("shutdown /a")
## [1] 1116
# Set options for analysis
use_mi = FALSE
shutdown = FALSE
report_ordinal = FALSE
do_priorsense = FALSE
get_bayesfactor = TRUE
check_models = TRUE #

if (get_bayesfactor) {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'BF', 'Rhat', 'ESS')
} else {
  stats_to_report <- c('CI', 'SE', 'pd', 'ROPE', 'Rhat', 'ESS')
}

options(
  dplyr.print_max = 100, 
  brms.backend = 'cmdstan',
  brms.file_refit = ifelse(use_mi, 'never', 'on_change'),
  brms.file_refit = 'on_change',
  #brms.file_refit = 'always',
  error = function() {
    beepr::beep(sound = 5)
    if (shutdown) {
      system("shutdown /s /t 180")
      quit(save = "no", status = 1)
    }
  }
  , es.use_symbols = TRUE
)


####################### Model parameters #######################

iterations = 12000 # 12'000 per chain to achieve 40'000
warmup = 2000 # 2000

# NO AR!!!
#corstr = 'ar'
#corstr = 'cosy_couple'
#corstr = 'cosy_couple:user'


################################################################

suffix_nopushing = paste0('_sensitivityNoPushing_', as.character(iterations))
suffix_onlypushing = paste0('_sensitivityOnlyPushing_', as.character(iterations))
df <- openxlsx::read.xlsx(file.path('long.xlsx'))
df_original <- df

df_double <- prepare_data(df, recode_pushing = FALSE, use_mi = use_mi)[[1]]

Constructing scales reshaping data (4field) centering data within and between

summary(df_double$pushing)

Min. 1st Qu. Median Mean 3rd Qu. Max. NA’s 0.0000 0.0000 0.0000 0.3468 0.0000 5.0000 2323

Modelling

# For indistinguishable Dyads
model_rows_fixed <- c(
    'Intercept', 
    # '-- WITHIN PERSON MAIN EFFECTS --', 
    'persuasion_self_cw', 
    'persuasion_partner_cw', 
    'pressure_self_cw', 
    'pressure_partner_cw', 
    'pushing_self_cw', 
    'pushing_partner_cw', 
    'day', 
    'plan_selfPlan',
    'plan_partnerPlan',
    'weartime_self_cw',
    
    # '-- BETWEEN PERSON MAIN EFFECTS',
    'persuasion_self_cb',
    'persuasion_partner_cb',
    'pressure_self_cb',
    'pressure_partner_cb',
    'pushing_self_cb',
    'pushing_partner_cb',
    'weartime_self_cb'
  )


model_rows_fixed_ordinal <- c(
  model_rows_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rows_fixed[2:length(model_rows_fixed)]
)

model_rows_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  'sd(persuasion_self_cw)',
  'sd(persuasion_partner_cw)',
  'sd(pressure_self_cw)',
  'sd(pressure_partner_cw)',
  'sd(pushing_self_cw)',
  'sd(pushing_partner_cw)',
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rows_random_ordinal <- c(model_rows_random,'disc')
# For indistinguishable Dyads
model_rownames_fixed <- c(
    "Intercept", 
    # "-- WITHIN PERSON MAIN EFFECTS --", 
    "Daily individual's experienced persuasion",  
    "Daily partner's experienced persuasion", 
    "Daily individual's experienced pressure", 
    "Daily partner's experienced pressure", 
    "Daily individual's experienced pushing", 
    "Daily partner's experienced pushing", 
    "Day", 
    "Own Actionplan",
    'Partner Actionplan',
    "Daily weartime",
    
    # "-- BETWEEN PERSON MAIN EFFECTS",
    "Mean individual's experienced persuasion", 
    "Mean partner's experienced persuasion", 
    "Mean individual's experienced pressure", 
    "Mean partner's experienced pressure", 
    "Mean individual's experienced pushing", 
    "Mean partner's experienced pushing", 
    "Mean weartime"
  )


model_rownames_fixed_ordinal <- c(
  model_rownames_fixed[1],
  'Intercept[1]',
  'Intercept[2]',
  'Intercept[3]',
  'Intercept[4]',
  'Intercept[5]',
  model_rownames_fixed[2:length(model_rownames_fixed)]
)

model_rownames_random <- c(
  # '--------------',
  # '-- RANDOM EFFECTS --',
  'sd(Intercept)', 
  "sd(Daily individual's experienced persuasion)", 
  "sd(Daily partner's experienced persuasion)", # OR partner received
  "sd(Daily individual's experienced pressure)", 
  "sd(Daily partner's experienced pressure)", 
  "sd(Daily individual's experienced pushing)", 
  "sd(Daily partner's experienced pushing)", 
  # '-- CORRELATION STRUCTURE -- ', 
  'sigma'
)

model_rownames_random_ordinal <- c(model_rownames_random,'disc')
rows_to_pack <- list(
  "Within-Person Effects" = c(2,11),
  "Between-Person Effects" = c(12,18),
  "Random Effects" = c(19, 25), 
  "Additional Parameters" = c(26,26)
  )


rows_to_pack_ordinal <- list(
  "Within-Person Effects" = c(2+5,11+5),
  "Between-Person Effects" = c(12+5,18+5),
  "Random Effects" = c(19+5, 25+5), 
  "Additional Parameters" = c(26+5,26+6)
  )

WITHOUT PUSHING

Self-Reported MVPA

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | dd | coupleID),
  
  hu = ~ persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | dd | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix_nopushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report, rope_range
## = rope_range_continuous, : Coefficients were exponentiated. Double check if
## this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 0.21*** 0.04 [ 0.15, 0.30] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 18379 22311 37.30*** 2.97 [31.85, 43.68] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 10134 18235
Within-Person Effects
Daily individual’s experienced persuasion 1.54*** 0.11 [ 1.35, 1.80] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 25035 23589 1.03 0.03 [ 0.98, 1.08] 0.841 [0.92, 1.08] 0.973 0.015 Very Strong Evidence for Null 1.000 16016 25099
Daily partner’s experienced persuasion 1.37*** 0.08 [ 1.23, 1.56] 1.000 [0.84, 1.20] 0.009 >100 Overwhelming Evidence 1.000 30330 26000 1.02 0.02 [ 0.98, 1.07] 0.835 [0.92, 1.08] 0.991 0.013 Very Strong Evidence for Null 1.000 18751 26921
Daily individual’s experienced pressure 1.02 0.15 [ 0.75, 1.39] 0.551 [0.84, 1.20] 0.765 0.075 Strong Evidence for Null 1.000 32143 26907 0.89* 0.05 [ 0.80, 0.99] 0.983 [0.92, 1.08] 0.256 0.080 Strong Evidence for Null 1.000 34259 33297
Daily partner’s experienced pressure 1.66** 0.34 [ 1.14, 2.76] 0.996 [0.84, 1.20] 0.042 3.558 Moderate Evidence 1.000 25056 18041 0.94 0.04 [ 0.85, 1.02] 0.928 [0.92, 1.08] 0.608 0.020 Very Strong Evidence for Null 1.000 33731 34026
Daily individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Day 0.80 0.12 [ 0.60, 1.07] 0.935 [0.84, 1.20] 0.378 0.236 Moderate Evidence for Null 1.000 47425 31069 0.96 0.06 [ 0.85, 1.08] 0.775 [0.92, 1.08] 0.681 0.009 Very Strong Evidence for Null 1.000 43576 31723
Own Actionplan 10.28*** 1.05 [ 8.41, 12.57] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 37997 31507 1.34*** 0.06 [ 1.22, 1.48] 1.000 [0.92, 1.08] 0.000 >100 Overwhelming Evidence 1.000 46497 30404
Partner Actionplan 1.23* 0.12 [ 1.01, 1.50] 0.982 [0.84, 1.20] 0.380 0.469 Weak Evidence for Null 1.000 39228 32059 1.07 0.05 [ 0.98, 1.16] 0.942 [0.92, 1.08] 0.600 0.025 Very Strong Evidence for Null 1.000 43919 30740
Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.65 0.57 [ 0.81, 3.29] 0.920 [0.84, 1.20] 0.151 0.490 Weak Evidence for Null 1.001 11799 19556 1.06 0.15 [ 0.80, 1.41] 0.666 [0.92, 1.08] 0.386 0.016 Very Strong Evidence for Null 1.000 10079 18377
Mean partner’s experienced persuasion 1.68 0.58 [ 0.83, 3.34] 0.926 [0.84, 1.20] 0.144 0.529 Weak Evidence for Null 1.000 11787 19210 1.04 0.15 [ 0.78, 1.37] 0.602 [0.92, 1.08] 0.413 0.015 Very Strong Evidence for Null 1.000 10081 17165
Mean individual’s experienced pressure 0.16*** 0.08 [ 0.06, 0.41] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 15084 23719 1.02 0.32 [ 0.55, 1.91] 0.527 [0.92, 1.08] 0.200 0.021 Very Strong Evidence for Null 1.001 6324 12200
Mean partner’s experienced pressure 0.34* 0.17 [ 0.12, 0.86] 0.989 [0.84, 1.20] 0.025 3.324 Moderate Evidence 1.000 15223 22711 0.84 0.26 [ 0.45, 1.56] 0.720 [0.92, 1.08] 0.167 0.024 Very Strong Evidence for Null 1.001 6356 11849
Mean individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.80 0.11 [0.61, 1.07] NA NA NA NA NA 1.000 15225 22010 0.31 0.04 [0.23, 0.41] NA NA NA NA NA 1.000 10568 17077
sd(Daily individual’s experienced persuasion) 0.24 0.09 [0.05, 0.43] NA NA NA NA NA 1.000 12645 10212 0.11 0.02 [0.07, 0.16] NA NA NA NA NA 1.000 17986 25110
sd(Daily partner’s experienced persuasion) 0.15 0.10 [0.01, 0.34] NA NA NA NA NA 1.000 10328 13402 0.08 0.02 [0.04, 0.13] NA NA NA NA NA 1.000 20247 19938
sd(Daily individual’s experienced pressure) 0.19 0.17 [0.01, 0.73] NA NA NA NA NA 1.000 16400 17101 0.08 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 13679 15950
sd(Daily partner’s experienced pressure) 0.34 0.27 [0.02, 1.14] NA NA NA NA NA 1.000 13168 17556 0.07 0.05 [0.00, 0.20] NA NA NA NA NA 1.000 16514 16084
sd(Daily individual’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.68 0.01 [0.66, 0.70] NA NA NA NA NA 1.000 43447 30427 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_hu_persuasion_partner_cb and
##   b_hu_persuasion_self_cb (r = 0.77), b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.79), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.85). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Device Based MVPA

Lognormal Model

formula <- bf(
  pa_obj ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +
    
    plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix_nopushing))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 109.79*** 5.82 [98.70, 122.10] 1.000 [0.94, 1.07] 0.000 >100 Overwhelming Evidence 1.001 8361 16470
Within-Person Effects
Daily individual’s experienced persuasion 1.03 0.01 [ 1.00, 1.06] 0.972 [0.94, 1.07] 0.992 0.034 Strong Evidence for Null 1.000 36290 31818
Daily partner’s experienced persuasion 1.02 0.02 [ 0.99, 1.05] 0.852 [0.94, 1.07] 0.999 0.010 Very Strong Evidence for Null 1.000 40350 33166
Daily individual’s experienced pressure 0.95 0.03 [ 0.89, 1.02] 0.937 [0.94, 1.07] 0.647 0.016 Very Strong Evidence for Null 1.000 60366 33325
Daily partner’s experienced pressure 0.98 0.03 [ 0.92, 1.05] 0.683 [0.94, 1.07] 0.929 0.005 Very Strong Evidence for Null 1.000 67097 32056
Daily individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Day 0.97 0.03 [ 0.91, 1.04] 0.770 [0.94, 1.07] 0.851 0.005 Very Strong Evidence for Null 1.000 91656 29083
Own Actionplan 1.07** 0.03 [ 1.02, 1.13] 0.997 [0.94, 1.07] 0.423 0.218 Moderate Evidence for Null 1.001 96483 27890
Partner Actionplan 1.06* 0.03 [ 1.01, 1.11] 0.987 [0.94, 1.07] 0.663 0.049 Strong Evidence for Null 1.000 93619 28077
Daily weartime 1.00*** 0.00 [ 1.00, 1.00] 1.000 [0.94, 1.07] 1.000 16.029 Strong Evidence 1.000 93707 29500
Between-Person Effects
Mean individual’s experienced persuasion 1.10 0.13 [ 0.86, 1.41] 0.791 [0.94, 1.07] 0.299 0.018 Very Strong Evidence for Null 1.000 8327 15077
Mean partner’s experienced persuasion 1.03 0.13 [ 0.81, 1.33] 0.604 [0.94, 1.07] 0.385 0.013 Very Strong Evidence for Null 1.000 8665 15408
Mean individual’s experienced pressure 0.90 0.12 [ 0.68, 1.18] 0.783 [0.94, 1.07] 0.269 0.011 Very Strong Evidence for Null 1.001 11838 20613
Mean partner’s experienced pressure 1.03 0.14 [ 0.79, 1.35] 0.588 [0.94, 1.07] 0.359 0.009 Very Strong Evidence for Null 1.000 10756 19400
Mean individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.928 [0.94, 1.07] 1.000 0.020 Very Strong Evidence for Null 1.000 48268 32570
Random Effects
sd(Intercept) 0.29 0.04 [0.23, 0.38] NA NA NA NA NA 1.000 11175 19228
sd(Daily individual’s experienced persuasion) 0.05 0.01 [0.02, 0.08] NA NA NA NA NA 1.000 23157 19605
sd(Daily partner’s experienced persuasion) 0.05 0.02 [0.03, 0.09] NA NA NA NA NA 1.000 20915 20170
sd(Daily individual’s experienced pressure) 0.04 0.04 [0.00, 0.15] NA NA NA NA NA 1.000 19919 24024
sd(Daily partner’s experienced pressure) 0.03 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 26364 21765
sd(Daily individual’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.58 0.01 [0.56, 0.59] NA NA NA NA NA 1.000 87822 27945 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_persuasion_partner_cb and
##   b_persuasion_self_cb (r = 0.87), b_pressure_self_cb and
##   b_persuasion_self_cb (r = 0.84), b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.77), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.76), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.86). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

# Nothing significant, no plots

Affect

Gaussian

formula <- bf(
  aff ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix_nopushing))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Warning in compute_rope(range = rope_range): Collinearity detected. Some VIFs
## are > 10. This may invalidate ROPE inferences!

Check for Multicollinearity

Low Correlation

              Term  VIF     VIF 95% CI Increased SE Tolerance
persuasion_self_cw 1.09 [ 1.06,  1.13]         1.04      0.92

persuasion_partner_cw 1.06 [ 1.04, 1.11] 1.03 0.94 pressure_self_cw 1.05 [ 1.03, 1.10] 1.03 0.95 pressure_partner_cw 1.04 [ 1.02, 1.09] 1.02 0.96 plan_self 1.29 [ 1.24, 1.34] 1.13 0.78 plan_partner 1.29 [ 1.24, 1.34] 1.13 0.78 day 1.04 [ 1.02, 1.09] 1.02 0.96 Tolerance 95% CI [0.89, 0.95] [0.90, 0.96] [0.91, 0.97] [0.92, 0.99] [0.75, 0.80] [0.75, 0.80] [0.91, 0.98]

Moderate Correlation

            Term  VIF     VIF 95% CI Increased SE Tolerance
pressure_self_cb 5.71 [ 5.40,  6.03]         2.39      0.18

pressure_partner_cb 5.54 [ 5.24, 5.85] 2.35 0.18 Tolerance 95% CI [0.17, 0.19] [0.17, 0.19]

High Correlation

              Term   VIF     VIF 95% CI Increased SE Tolerance
persuasion_self_cb 12.17 [11.48, 12.90]         3.49      0.08

persuasion_partner_cb 12.05 [11.37, 12.77] 3.47 0.08 Tolerance 95% CI [0.08, 0.09] [0.08, 0.09]

## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.65*** 0.11 [ 3.44, 3.86] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.001 4597 9596
Within-Person Effects
Daily individual’s experienced persuasion 0.00 0.02 [-0.04, 0.05] 0.545 [-0.11, 0.11] 1.000 0.004 Very Strong Evidence for Null 1.000 43828 29691
Daily partner’s experienced persuasion 0.03 0.02 [-0.01, 0.08] 0.909 [-0.11, 0.11] 0.999 0.010 Very Strong Evidence for Null 1.000 28834 30041
Daily individual’s experienced pressure -0.02 0.05 [-0.13, 0.08] 0.679 [-0.11, 0.11] 0.952 0.004 Very Strong Evidence for Null 1.000 52388 29167
Daily partner’s experienced pressure 0.02 0.05 [-0.09, 0.12] 0.634 [-0.11, 0.11] 0.959 0.004 Very Strong Evidence for Null 1.000 49735 28360
Daily individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Day 0.25*** 0.06 [ 0.14, 0.36] 1.000 [-0.11, 0.11] 0.008 24.739 Strong Evidence 1.000 70186 29727
Own Actionplan 0.12** 0.04 [ 0.04, 0.20] 0.999 [-0.11, 0.11] 0.448 0.338 Weak Evidence for Null 1.000 82577 28323
Partner Actionplan -0.01 0.04 [-0.09, 0.06] 0.621 [-0.11, 0.11] 0.995 0.003 Very Strong Evidence for Null 1.000 73169 27432
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.42 0.26 [-0.10, 0.93] 0.947 [-0.11, 0.11] 0.098 0.051 Strong Evidence for Null 1.001 4951 10892
Mean partner’s experienced persuasion 0.33 0.26 [-0.19, 0.84] 0.895 [-0.11, 0.11] 0.162 0.030 Very Strong Evidence for Null 1.001 4940 10928
Mean individual’s experienced pressure -0.35 0.28 [-0.90, 0.20] 0.897 [-0.11, 0.11] 0.147 0.021 Very Strong Evidence for Null 1.001 5879 12671
Mean partner’s experienced pressure -0.28 0.28 [-0.83, 0.27] 0.844 [-0.11, 0.11] 0.194 0.016 Very Strong Evidence for Null 1.001 5842 12994
Mean individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.60 0.07 [0.48, 0.78] NA NA NA NA NA 1.000 8044 15302
sd(Daily individual’s experienced persuasion) 0.04 0.03 [0.00, 0.11] NA NA NA NA NA 1.000 10880 14460
sd(Daily partner’s experienced persuasion) 0.08 0.03 [0.02, 0.14] NA NA NA NA NA 1.000 11899 8164
sd(Daily individual’s experienced pressure) 0.06 0.06 [0.00, 0.23] NA NA NA NA NA 1.000 17936 21028
sd(Daily partner’s experienced pressure) 0.07 0.06 [0.00, 0.24] NA NA NA NA NA 1.000 19481 20163
sd(Daily individual’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma 0.96 0.01 [0.94, 0.98] NA NA NA NA NA 1.000 69129 30080 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()
## Possible multicollinearity between b_pressure_partner_cb and
##   b_persuasion_self_cb (r = 0.86), b_pressure_self_cb and
##   b_persuasion_partner_cb (r = 0.86), b_pressure_partner_cb and
##   b_pressure_self_cb (r = 0.82). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Reactance

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix_nopushing))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 3.37*** 0.98 [ 1.93, 6.02] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 35331 31342
Intercept[2] 7.16*** 2.15 [ 4.03, 13.08] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 34683 31744
Intercept[3] 18.65*** 6.00 [ 10.11, 35.29] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 35222 31270
Intercept[4] 77.42*** 28.66 [ 38.09, 163.46] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 39555 32044
Intercept[5] 2234.28*** 1447.67 [696.10, 9042.07] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 49254 29590
Within-Person Effects
Daily individual’s experienced persuasion 0.85* 0.07 [ 0.72, 1.00] 0.975 [0.84, 1.20] 0.598 0.334 Weak Evidence for Null 1.000 37527 29881
Daily partner’s experienced persuasion 1.01 0.09 [ 0.83, 1.21] 0.548 [0.84, 1.20] 0.942 0.044 Strong Evidence for Null 1.000 37362 29459
Daily individual’s experienced pressure 2.01** 0.34 [ 1.34, 2.75] 0.999 [0.84, 1.20] 0.008 8.308 Moderate Evidence 1.000 26737 24583
Daily partner’s experienced pressure 1.12 0.23 [ 0.70, 1.79] 0.710 [0.84, 1.20] 0.537 0.055 Strong Evidence for Null 1.000 29868 23638
Daily individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Day 1.74 0.59 [ 0.89, 3.40] 0.947 [0.84, 1.20] 0.120 0.136 Moderate Evidence for Null 1.000 50402 30696
Own Actionplan 1.07 0.25 [ 0.67, 1.72] 0.606 [0.84, 1.20] 0.534 0.037 Strong Evidence for Null 1.000 53488 29621
Partner Actionplan 0.87 0.20 [ 0.55, 1.36] 0.733 [0.84, 1.20] 0.483 0.043 Strong Evidence for Null 1.000 55337 30868
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 0.96 0.48 [ 0.36, 2.62] 0.530 [0.84, 1.20] 0.284 0.056 Strong Evidence for Null 1.000 19529 26959
Mean partner’s experienced persuasion 0.77 0.41 [ 0.27, 2.30] 0.689 [0.84, 1.20] 0.231 0.062 Strong Evidence for Null 1.000 21295 26278
Mean individual’s experienced pressure 3.52* 1.92 [ 1.20, 10.59] 0.989 [0.84, 1.20] 0.021 0.750 Weak Evidence for Null 1.000 22874 28186
Mean partner’s experienced pressure 1.27 0.73 [ 0.39, 3.93] 0.659 [0.84, 1.20] 0.224 0.055 Strong Evidence for Null 1.000 20597 27389
Mean individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.78 0.19 [0.45, 1.22] NA NA NA NA NA 1.000 15053 24788
sd(Daily individual’s experienced persuasion) 0.19 0.13 [0.01, 0.46] NA NA NA NA NA 1.000 7448 14855
sd(Daily partner’s experienced persuasion) 0.19 0.13 [0.01, 0.48] NA NA NA NA NA 1.000 13143 16496
sd(Daily individual’s experienced pressure) 0.42 0.23 [0.04, 0.97] NA NA NA NA NA 1.000 10254 12322
sd(Daily partner’s experienced pressure) 0.30 0.30 [0.01, 1.38] NA NA NA NA NA 1.001 10915 20948
sd(Daily individual’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.8), b_Intercept[4] and b_Intercept[3] (r = 0.85), b_pressure_self_cb
##   and b_persuasion_self_cb (r = 0.81), b_pressure_partner_cb and
##   b_persuasion_partner_cb (r = 0.85). This might lead to inappropriate
##   results. See 'Details' in '?rope'.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw +

    persuasion_self_cb + persuasion_partner_cb +
    pressure_self_cb + pressure_partner_cb +

    plan_self + plan_partner +
    day +
    
    # Random effects
    (persuasion_self_cw + persuasion_partner_cw +
    pressure_self_cw + pressure_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix_nopushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
summary_is_reactance <- summarize_brms(
  is_reactance,
  stats_to_report = stats_to_report, 
  model_rows_fixed = model_rows_fixed, 
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.34** 0.11 [0.18, 0.64] 0.999 [0.83, 1.20] 0.003 11.618 Strong Evidence 1.000 35375 31276
Within-Person Effects
Daily individual’s experienced persuasion 0.86 0.08 [0.70, 1.03] 0.952 [0.83, 1.20] 0.618 0.209 Moderate Evidence for Null 1.000 36196 32792
Daily partner’s experienced persuasion 1.12 0.15 [0.86, 1.50] 0.799 [0.83, 1.20] 0.687 0.085 Strong Evidence for Null 1.000 26249 29181
Daily individual’s experienced pressure 2.09* 0.55 [1.22, 4.32] 0.993 [0.83, 1.20] 0.019 2.746 Weak Evidence 1.000 24026 20657
Daily partner’s experienced pressure 1.30 0.53 [0.55, 4.09] 0.756 [0.83, 1.20] 0.286 0.115 Moderate Evidence for Null 1.000 20941 19405
Daily individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Day 1.94 0.73 [0.92, 4.08] 0.959 [0.83, 1.20] 0.090 0.190 Moderate Evidence for Null 1.000 63520 30318
Own Actionplan 1.14 0.33 [0.65, 2.04] 0.679 [0.83, 1.20] 0.424 0.049 Strong Evidence for Null 1.000 59668 30821
Partner Actionplan 0.85 0.24 [0.49, 1.47] 0.720 [0.83, 1.20] 0.419 0.052 Strong Evidence for Null 1.000 66858 29417
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.65 1.06 [0.47, 6.38] 0.788 [0.83, 1.20] 0.166 0.097 Strong Evidence for Null 1.000 15884 25354
Mean partner’s experienced persuasion 1.08 0.74 [0.28, 4.61] 0.543 [0.83, 1.20] 0.207 0.071 Strong Evidence for Null 1.000 15991 24607
Mean individual’s experienced pressure 25.87*** 27.32 [3.74, 268.79] 1.000 [0.83, 1.20] 0.001 28.786 Strong Evidence 1.000 18246 24515
Mean partner’s experienced pressure 0.86 0.98 [0.08, 7.19] 0.553 [0.83, 1.20] 0.126 0.102 Moderate Evidence for Null 1.000 14827 23586
Mean individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.27 0.26 [0.84, 1.90] NA NA NA NA NA 1.000 13250 22431
sd(Daily individual’s experienced persuasion) 0.22 0.15 [0.02, 0.53] NA NA NA NA NA 1.000 7931 13476
sd(Daily partner’s experienced persuasion) 0.44 0.17 [0.14, 0.87] NA NA NA NA NA 1.001 15843 14382
sd(Daily individual’s experienced pressure) 0.71 0.54 [0.04, 2.07] NA NA NA NA NA 1.000 6881 13518
sd(Daily partner’s experienced pressure) 0.78 0.68 [0.04, 2.79] NA NA NA NA NA 1.000 11946 19336
sd(Daily individual’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pushing) NA NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix_nopushing, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)
## 
## Attaching package: 'rvest'
## The following object is masked from 'package:readr':
## 
##     guess_encoding
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 0.21*** [ 0.15, 0.30] 37.30*** [31.85, 43.68] 109.79*** [98.70, 122.10] 3.65*** [ 3.44, 3.86] 0.34** [0.18, 0.64]
Within-Person Effects
Daily individual’s experienced persuasion 1.54*** [ 1.35, 1.80] 1.03 [ 0.98, 1.08] 1.03 [ 1.00, 1.06] 0.00 [-0.04, 0.05] 0.86 [0.70, 1.03]
Daily partner’s experienced persuasion 1.37*** [ 1.23, 1.56] 1.02 [ 0.98, 1.07] 1.02 [ 0.99, 1.05] 0.03 [-0.01, 0.08] 1.12 [0.86, 1.50]
Daily individual’s experienced pressure 1.02 [ 0.75, 1.39] 0.89* [ 0.80, 0.99] 0.95 [ 0.89, 1.02] -0.02 [-0.13, 0.08] 2.09* [1.22, 4.32]
Daily partner’s experienced pressure 1.66** [ 1.14, 2.76] 0.94 [ 0.85, 1.02] 0.98 [ 0.92, 1.05] 0.02 [-0.09, 0.12] 1.30 [0.55, 4.09]
Daily individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA
Day 0.80 [ 0.60, 1.07] 0.96 [ 0.85, 1.08] 0.97 [ 0.91, 1.04] 0.25*** [ 0.14, 0.36] 1.94 [0.92, 4.08]
Own Actionplan 10.28*** [ 8.41, 12.57] 1.34*** [ 1.22, 1.48] 1.07** [ 1.02, 1.13] 0.12** [ 0.04, 0.20] 1.14 [0.65, 2.04]
Partner Actionplan 1.23* [ 1.01, 1.50] 1.07 [ 0.98, 1.16] 1.06* [ 1.01, 1.11] -0.01 [-0.09, 0.06] 0.85 [0.49, 1.47]
Daily weartime NA NA NA NA 1.00*** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion 1.65 [ 0.81, 3.29] 1.06 [ 0.80, 1.41] 1.10 [ 0.86, 1.41] 0.42 [-0.10, 0.93] 1.65 [0.47, 6.38]
Mean partner’s experienced persuasion 1.68 [ 0.83, 3.34] 1.04 [ 0.78, 1.37] 1.03 [ 0.81, 1.33] 0.33 [-0.19, 0.84] 1.08 [0.28, 4.61]
Mean individual’s experienced pressure 0.16*** [ 0.06, 0.41] 1.02 [ 0.55, 1.91] 0.90 [ 0.68, 1.18] -0.35 [-0.90, 0.20] 25.87*** [3.74, 268.79]
Mean partner’s experienced pressure 0.34* [ 0.12, 0.86] 0.84 [ 0.45, 1.56] 1.03 [ 0.79, 1.35] -0.28 [-0.83, 0.27] 0.86 [0.08, 7.19]
Mean individual’s experienced pushing NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pushing NA NA NA NA NA NA NA NA NA NA
Mean weartime NA NA NA NA 1.00 [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 0.80 [0.61, 1.07] 0.31 [0.23, 0.41] 0.29 [0.23, 0.38] 0.60 [0.48, 0.78] 1.27 [0.84, 1.90]
sd(Daily individual’s experienced persuasion) 0.24 [0.05, 0.43] 0.11 [0.07, 0.16] 0.05 [0.02, 0.08] 0.04 [0.00, 0.11] 0.22 [0.02, 0.53]
sd(Daily partner’s experienced persuasion) 0.15 [0.01, 0.34] 0.08 [0.04, 0.13] 0.05 [0.03, 0.09] 0.08 [0.02, 0.14] 0.44 [0.14, 0.87]
sd(Daily individual’s experienced pressure) 0.19 [0.01, 0.73] 0.08 [0.00, 0.24] 0.04 [0.00, 0.15] 0.06 [0.00, 0.23] 0.71 [0.04, 2.07]
sd(Daily partner’s experienced pressure) 0.34 [0.02, 1.14] 0.07 [0.00, 0.20] 0.03 [0.00, 0.11] 0.07 [0.00, 0.24] 0.78 [0.04, 2.79]
sd(Daily individual’s experienced pushing) NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pushing) NA NA NA NA NA NA NA NA NA NA
Additional Parameters
sigma NA NA 0.68 [0.66, 0.70] 0.58 [0.56, 0.59] 0.96 [0.94, 0.98] NA NA

ONLY PUSHING

Hurdle Lognormal Model

formula <- bf(
  pa_sub ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID),
  
  hu = ~ pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day + 
    
    # Random effects
    (1 + pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
) 

prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "b", dpar = "hu")
  , brms::set_prior("normal(0, 50)", class = "Intercept") # for non-zero PA
  , brms::set_prior("normal(0.5, 2.5)", class = "Intercept", dpar = 'hu') # hurdle part
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = hurdle_lognormal()
#)

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_sub <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::hurdle_lognormal(), 
  #family = brms::hurdle_negbinomial(), 
  #family = brms::hurdle_poisson(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 42,
  file = file.path("models_cache_brms", paste0("pa_sub_hu_lognormal", suffix_onlypushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
pa_sub_digest <- digest::digest(pa_sub)
# rope range for continuous part of the model
rope_factor <- sd(log(pa_sub$data$pa_sub[pa_sub$data$pa_sub > 0]))
rope_range_continuous = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_sub <- summarize_brms(
  pa_sub, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_continuous,
  hu_rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_sub, stats_to_report = stats_to_report, rope_range
## = rope_range_continuous, : Coefficients were exponentiated. Double check if
## this was intended.
# Print the updated dataframe
summary_pa_sub %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.)_hu SE_hu 95% CI_hu pd_hu ROPE_hu inside ROPE_hu BF_hu BF_Evidence_hu Rhat_hu Bulk_ESS_hu Tail_ESS_hu exp(Est.)_nonzero SE_nonzero 95% CI_nonzero pd_nonzero ROPE_nonzero inside ROPE_nonzero BF_nonzero BF_Evidence_nonzero Rhat_nonzero Bulk_ESS_nonzero Tail_ESS_nonzero
Intercept 3.79*** 0.97 [ 2.28, 6.37] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.001 9084 17256 57.44*** 4.50 [49.33, 67.17] 1.000 [0.93, 1.08] 0.000 >100 Overwhelming Evidence 1.000 11116 19378
Within-Person Effects
Daily individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pushing 1.55** 0.23 [ 1.17, 2.18] 0.998 [0.84, 1.20] 0.033 6.417 Moderate Evidence 1.000 24317 24788 0.97 0.03 [ 0.91, 1.03] 0.852 [0.93, 1.08] 0.934 0.017 Very Strong Evidence for Null 1.000 33376 26045
Daily partner’s experienced pushing 1.84*** 0.30 [ 1.36, 2.72] 1.000 [0.84, 1.20] 0.004 54.092 Very Strong Evidence 1.000 21864 23974 0.97 0.03 [ 0.90, 1.03] 0.851 [0.93, 1.08] 0.904 0.018 Very Strong Evidence for Null 1.000 27989 26716
Day 0.61* 0.15 [ 0.37, 0.99] 0.978 [0.84, 1.20] 0.100 0.892 Weak Evidence for Null 1.000 59962 30420 0.91 0.07 [ 0.79, 1.05] 0.896 [0.93, 1.08] 0.399 0.020 Very Strong Evidence for Null 1.000 66158 30687
Own Actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pushing 0.47** 0.13 [ 0.25, 0.80] 0.998 [0.84, 1.20] 0.017 7.148 Moderate Evidence 1.000 19631 23385 1.32* 0.18 [ 1.01, 1.72] 0.977 [0.93, 1.08] 0.068 0.165 Moderate Evidence for Null 1.000 9026 16236
Mean partner’s experienced pushing 0.63 0.18 [ 0.34, 1.08] 0.955 [0.84, 1.20] 0.142 0.568 Weak Evidence for Null 1.000 20071 23775 1.17 0.16 [ 0.89, 1.53] 0.881 [0.93, 1.08] 0.228 0.041 Strong Evidence for Null 1.000 8988 15578
Mean weartime NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.26 0.20 [0.92, 1.75] NA NA NA NA NA 1.000 11164 19188 0.37 0.05 [0.28, 0.50] NA NA NA NA NA 1.000 10863 20431
sd(Daily individual’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pushing) 0.47 0.24 [0.07, 1.02] NA NA NA NA NA 1.000 7788 9454 0.06 0.04 [0.00, 0.16] NA NA NA NA NA 1.000 7540 11307
sd(Daily partner’s experienced pushing) 0.51 0.24 [0.09, 1.10] NA NA NA NA NA 1.000 8390 8477 0.09 0.04 [0.01, 0.18] NA NA NA NA NA 1.000 7899 9193
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 0.67 0.02 [0.64, 0.70] NA NA NA NA NA 1.000 41396 29445 
# Plot continuous part of model

variable <- c(
  '(Intercept)',
  'b_persuasion_self_cw',
  'b_persuasion_partner_cw',
  'b_pressure_self_cw',
  'b_pressure_partner_cw',
  'b_pushing_self_cw',
  'b_pushing_partner_cw'
)


plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

plot(
  bayestestR::rope(
    pa_sub, 
    parameter = variable, 
    range = rope_range_continuous,
    verbose = F,
    ci = 1
  )
) + theme_bw()

# Hurdle part of the model
variable <- c(
  'b_hu_persuasion_self_cw',
  'b_hu_persuasion_partner_cw',
  'b_hu_pressure_self_cw',
  'b_hu_pressure_partner_cw',
  'b_hu_pushing_self_cw',
  'b_hu_pushing_partner_cw'
)

plot(
  bayestestR::p_direction(pa_sub, parameter = variable),
  priors = TRUE
) + theme_bw()

# The rope range for the bernoulli part of the model is -0.18, 0.18
plot(
  bayestestR::rope(pa_sub, parameter = variable, range = c(-0.18, 0.18), ci = 1),
  verbose = FALSE
) + theme_bw()
## Possible multicollinearity between b_pushing_partner_cb and
##   b_pushing_self_cb (r = 0.77). This might lead to inappropriate results.
##   See 'Details' in '?rope'.

Device Based MVPA

Lognormal Model

formula <- bf(
  pa_obj ~ 
   
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day + weartime_self_cw + weartime_self_cb +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 50)", class = "Intercept") 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = lognormal()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

pa_obj_log <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = lognormal(),
  #control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("pa_obj_log_gaussian", suffix_onlypushing))
)
## Warning: Rows containing NAs were excluded from the model.
pa_obj_log_digest <- digest::digest(pa_obj_log)
# rope range for lognormal model
rope_factor <- sd(log(pa_obj_log$data$pa_obj))
rope_range_log = c(-0.1 * rope_factor, 0.1 * rope_factor)

summary_pa_obj <- summarize_brms(
  pa_obj_log, 
  stats_to_report = stats_to_report,
  rope_range = rope_range_log,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
## Warning in summarize_brms(pa_obj_log, stats_to_report = stats_to_report, :
## Coefficients were exponentiated. Double check if this was intended.
summary_pa_obj %>%
  print_df(rows_to_pack = rows_to_pack)
exp(Est.) SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 124.09*** 7.54 [109.75, 140.55] 1.000 [0.94, 1.06] 0.000 >100 Overwhelming Evidence 1.001 6735 13275
Within-Person Effects
Daily individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pushing 1.02 0.03 [ 0.97, 1.08] 0.806 [0.94, 1.06] 0.910 0.014 Very Strong Evidence for Null 1.000 19888 22821
Daily partner’s experienced pushing 1.02 0.02 [ 0.97, 1.06] 0.789 [0.94, 1.06] 0.978 0.010 Very Strong Evidence for Null 1.000 35178 27905
Day 0.97 0.05 [ 0.87, 1.08] 0.717 [0.94, 1.06] 0.662 0.008 Very Strong Evidence for Null 1.000 45016 30794
Own Actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily weartime 1.00** 0.00 [ 1.00, 1.00] 0.998 [0.94, 1.06] 1.000 0.377 Weak Evidence for Null 1.000 47851 29455
Between-Person Effects
Mean individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pushing 0.97 0.06 [ 0.85, 1.10] 0.688 [0.94, 1.06] 0.608 0.011 Very Strong Evidence for Null 1.000 15998 25391
Mean partner’s experienced pushing 1.04 0.07 [ 0.92, 1.18] 0.736 [0.94, 1.06] 0.579 0.012 Very Strong Evidence for Null 1.000 13600 22072
Mean weartime 1.00 0.00 [ 1.00, 1.00] 0.813 [0.94, 1.06] 1.000 0.016 Very Strong Evidence for Null 1.000 18544 26401
Random Effects
sd(Intercept) 0.31 0.04 [0.24, 0.42] NA NA NA NA NA 1.000 8512 14572
sd(Daily individual’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pushing) 0.09 0.04 [0.02, 0.17] NA NA NA NA NA 1.001 8388 7348
sd(Daily partner’s experienced pushing) 0.03 0.03 [0.00, 0.10] NA NA NA NA NA 1.001 12784 16579
Additional Parameters
sigma 0.54 0.01 [0.52, 0.57] NA NA NA NA NA 1.000 37252 30307 
plot(
  bayestestR::p_direction(pa_obj_log),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(pa_obj_log, range = rope_range_log, ci = 1)
) + theme_bw()

# Nothing significant, no plots

Affect

Gaussian

formula <- bf(
  aff ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 5)", class = "b")
  ,brms::set_prior("normal(0, 20)", class = "Intercept", lb=1, ub=6)
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
  , brms::set_prior("student_t(3, 0, 2.5)", class = "sigma", lb = 0)
)

#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = gaussian()
#  )

#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

mood_gauss <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = gaussian(),
  #control = list(adapt_delta = 0.95, max_treedepth = 15),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("mood_gauss", suffix_onlypushing))
)
## Warning: Rows containing NAs were excluded from the model.
mood_gauss_digest <- digest::digest(mood_gauss)
summary_mood <- summarize_brms(
  mood_gauss, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = F) 
## Sampling priors, please wait...
summary_mood %>%
  print_df(rows_to_pack = rows_to_pack)
Est. SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 3.77*** 0.11 [ 3.54, 4.00] 1.000 [-0.11, 0.11] 0.000 >100 Overwhelming Evidence 1.000 4876 8494
Within-Person Effects
Daily individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pushing 0.02 0.03 [-0.05, 0.09] 0.691 [-0.11, 0.11] 0.994 0.006 Very Strong Evidence for Null 1.000 28231 27942
Daily partner’s experienced pushing 0.09* 0.04 [ 0.00, 0.17] 0.978 [-0.11, 0.11] 0.700 0.063 Strong Evidence for Null 1.000 15128 22181
Day 0.28** 0.09 [ 0.11, 0.45] 0.999 [-0.11, 0.11] 0.024 0.803 Weak Evidence for Null 1.000 43377 32307
Own Actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pushing 0.13 0.10 [-0.08, 0.33] 0.894 [-0.11, 0.11] 0.414 0.017 Very Strong Evidence for Null 1.000 14612 23364
Mean partner’s experienced pushing 0.09 0.10 [-0.11, 0.29] 0.817 [-0.11, 0.11] 0.543 0.012 Very Strong Evidence for Null 1.000 13156 21082
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.61 0.08 [0.48, 0.80] NA NA NA NA NA 1.001 6607 13959
sd(Daily individual’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pushing) 0.07 0.04 [0.01, 0.16] NA NA NA NA NA 1.000 13468 11788
sd(Daily partner’s experienced pushing) 0.14 0.05 [0.06, 0.24] NA NA NA NA NA 1.000 15028 11785
Additional Parameters
sigma 0.91 0.02 [0.87, 0.94] NA NA NA NA NA 1.000 44921 30740 
plot(
  bayestestR::p_direction(mood_gauss),
  priors = TRUE
)  + 
  coord_cartesian(xlim = c(-3, 3)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(mood_gauss, ci = 1)
) + theme_bw()

conditional_spaghetti(
  mood_gauss, 
  effects = c('pushing_partner_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_partner_cw

Reactance

Ordinal

df_double$reactance_ordinal <- factor(df_double$reactance,
                                      levels = 0:5, 
                                      ordered = TRUE)

formula <- bf(
  reactance_ordinal ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
)


prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = cumulative() # HURDLE_CUMULATIVE
#  )


#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

reactance_ordinal <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::cumulative(),
  control = list(adapt_delta = 0.90),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777
  , file = file.path("models_cache_brms", paste0("reactance_ordinal", suffix_onlypushing))
)
## Warning: Rows containing NAs were excluded from the model.
reactance_ordinal_digest <- digest::digest(reactance_ordinal)
summary_reactance_ordinal <- summarize_brms(
  reactance_ordinal, 
  stats_to_report = stats_to_report,
  rope_range = c(-0.18, 0.18),
  model_rows_fixed = model_rows_fixed_ordinal,
  model_rows_random = model_rows_random_ordinal,
  model_rownames_fixed = model_rownames_fixed_ordinal,
  model_rownames_random = model_rownames_random_ordinal,
  exponentiate = T) 
## Sampling priors, please wait...
summary_reactance_ordinal %>%
  print_df(rows_to_pack = rows_to_pack_ordinal)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept NA NA NA NA NA NA NA NA NA NA NA
Intercept[1] 5.38*** 1.58 [ 3.07, 9.91] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 50877 32738
Intercept[2] 9.82*** 3.06 [ 5.44, 18.77] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 51473 32534
Intercept[3] 22.93*** 7.86 [ 11.99, 46.72] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 54044 33339
Intercept[4] 81.33*** 34.20 [ 37.44, 196.14] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 61913 31792
Intercept[5] 768.96*** 600.61 [201.43, 4928.66] 1.000 [0.84, 1.20] 0.000 >100 Overwhelming Evidence 1.000 79350 28894
Within-Person Effects
Daily individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pushing 1.30* 0.14 [ 1.04, 1.62] 0.988 [0.84, 1.20] 0.225 0.869 Weak Evidence for Null 1.000 61346 32371
Daily partner’s experienced pushing 0.97 0.13 [ 0.74, 1.25] 0.603 [0.84, 1.20] 0.812 0.057 Strong Evidence for Null 1.000 60592 28030
Day 1.62 0.76 [ 0.65, 4.05] 0.850 [0.84, 1.20] 0.181 0.092 Strong Evidence for Null 1.000 76425 31582
Own Actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pushing 2.30** 0.63 [ 1.36, 4.03] 0.999 [0.84, 1.20] 0.007 6.984 Moderate Evidence 1.000 55011 33534
Mean partner’s experienced pushing 1.24 0.35 [ 0.71, 2.15] 0.782 [0.84, 1.20] 0.364 0.075 Strong Evidence for Null 1.000 59497 33931
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 0.85 0.21 [0.48, 1.35] NA NA NA NA NA 1.000 21301 28454
sd(Daily individual’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pushing) 0.17 0.13 [0.01, 0.50] NA NA NA NA NA 1.000 15953 21118
sd(Daily partner’s experienced pushing) 0.16 0.14 [0.01, 0.58] NA NA NA NA NA 1.000 20500 24746
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA
disc 1.00 0.00 [1.00, 1.00] NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(reactance_ordinal),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()

plot(
  bayestestR::rope(reactance_ordinal, range = c(-0.18, 0.18), ci = 1)
) + theme_bw()
## Possible multicollinearity between b_Intercept[4] and b_Intercept[2] (r
##   = 0.75), b_Intercept[4] and b_Intercept[3] (r = 0.81). This might lead
##   to inappropriate results. See 'Details' in '?rope'.

Binary

introduce_binary_reactance <- function(data) {
  data$is_reactance <- factor(data$reactance > 0, levels = c(FALSE, TRUE), labels = c(0, 1))
  return(data)
}



df_double <- introduce_binary_reactance(df_double)
if (use_mi) {
  for (i in seq_along(implist)) {
    implist[[i]] <- introduce_binary_reactance(implist[[i]])
  }
}


formula <- bf(
  is_reactance ~ 
    pushing_self_cw + pushing_partner_cw +
    pushing_self_cb + pushing_partner_cb +
    
    #plan_self + plan_partner +
    day +
    
    # Random effects
    (pushing_self_cw + pushing_partner_cw | coupleID)
  
  , decomp = 'QR'
  #, autocor = autocor_str
  )



prior1 <- c(
  brms::set_prior("normal(0, 2.5)", class = "b")
  , brms::set_prior("normal(0, 10)", class = "Intercept", lb=0, ub=5) 
  , brms::set_prior("normal(0, 2)", class = "sd", group = "coupleID", lb = 0)
)


#brms::validate_prior(
#  prior1, 
#  formula = formula, 
#  data = df_double, 
#  family = bernoulli()
#  )



#df_minimal <- df_double[, c("AorB", all.vars(as.formula(formula)))]

is_reactance <- my_brm(
  mi = use_mi, 
  imputed_data = implist,
  
  formula = formula, 
  prior = prior1,
  data = df_double, 
  family = brms::bernoulli(),
  control = list(adapt_delta = 0.95),
  iter = iterations,
  warmup = warmup,
  chains = 4,
  cores = 4,
  seed = 7777,
  file = file.path("models_cache_brms", paste0("is_reactance", suffix_onlypushing))
  #, file_refit = 'always'
)
## Warning: Rows containing NAs were excluded from the model.
is_reactance_digest <- digest::digest(is_reactance)
summary_is_reactance <- summarize_brms(
  is_reactance, 
  stats_to_report = stats_to_report,
  model_rows_fixed = model_rows_fixed,
  model_rows_random = model_rows_random,
  model_rownames_fixed = model_rownames_fixed,
  model_rownames_random = model_rownames_random,
  exponentiate = T) 
## Sampling priors, please wait...
summary_is_reactance %>%
  print_df(rows_to_pack = rows_to_pack)
OR SE 95% CI pd ROPE inside ROPE BF BF_Evidence Rhat Bulk_ESS Tail_ESS
Intercept 0.40*** 0.11 [0.22, 0.67] 1.000 [0.83, 1.20] 0.003 19.245 Strong Evidence 1.000 28966 29532
Within-Person Effects
Daily individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pushing 1.36* 0.18 [1.06, 1.84] 0.992 [0.83, 1.20] 0.162 1.161 Weak Evidence 1.000 31218 25187
Daily partner’s experienced pushing 1.15 0.26 [0.79, 2.07] 0.752 [0.83, 1.20] 0.518 0.103 Moderate Evidence for Null 1.000 19898 22298
Day 1.78 0.90 [0.68, 4.81] 0.879 [0.83, 1.20] 0.152 0.114 Moderate Evidence for Null 1.000 58577 32613
Own Actionplan NA NA NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA NA NA NA NA NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pushing 5.68*** 2.76 [2.38, 16.59] 1.000 [0.83, 1.20] 0.000 >100 Overwhelming Evidence 1.000 16786 24262
Mean partner’s experienced pushing 2.41 1.23 [0.94, 7.27] 0.965 [0.83, 1.20] 0.062 0.465 Weak Evidence for Null 1.000 16636 22004
Mean weartime NA NA NA NA NA NA NA NA NA NA NA
Random Effects
sd(Intercept) 1.53 0.30 [1.04, 2.25] NA NA NA NA NA 1.000 16752 25783
sd(Daily individual’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pressure) NA NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pushing) 0.17 0.14 [0.01, 0.55] NA NA NA NA NA 1.000 13633 20559
sd(Daily partner’s experienced pushing) 0.41 0.28 [0.03, 1.12] NA NA NA NA NA 1.000 12144 15337
Additional Parameters
sigma NA NA NA NA NA NA NA NA NA NA NA 
plot(
  bayestestR::p_direction(is_reactance),
  priors = TRUE
) + 
  coord_cartesian(xlim = c(-6, 6)) +
  theme_bw()
## Warning in `==.default`(dens$Parameter, parameter): longer object length is not
## a multiple of shorter object length
## Warning in is.na(e1) | is.na(e2): longer object length is not a multiple of
## shorter object length

plot(
  bayestestR::rope(is_reactance, ci = 1)
) + theme_bw()

conditional_spaghetti(
  is_reactance, 
  effects = c('pushing_self_cw'),
  group_var = 'coupleID',
  plot_full_range = TRUE
)

$pushing_self_cw

Report All Models

summary_all_models <- report_side_by_side(
  pa_sub,
  pa_obj_log,
  mood_gauss,
  is_reactance,
  
  stats_to_report = c('CI'),
  model_rows_random = model_rows_random,
  model_rows_fixed = model_rows_fixed,
  model_rownames_random = model_rownames_random,
  model_rownames_fixed = model_rownames_fixed
)

[1] “pa_sub”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “pa_obj_log”

## Warning in summarize_brms(model, exponentiate = exponentiate, stats_to_report =
## stats_to_report, : Coefficients were exponentiated. Double check if this was
## intended.

[1] “mood_gauss” [1] “is_reactance”

summary_all_models <- summary_all_models %>%
  print_df(rows_to_pack = rows_to_pack) %>%
  add_header_above(
    c(
      " ", "Hurdle Component" = 2, "Non-Zero Component" = 2,
      " " = 6
    )
  ) %>%
  add_header_above(
    c(" ", "Subjective MVPA Hurdle Lognormal" = 4,  
      "Device-Based MVPA Log (Gaussian)" = 2, 
      "Mood Gaussian" = 2,
      #"Reactance Ordinal" = 2,
      "Reactance Dichotome" = 2
    )
  )

export_xlsx(
  summary_all_models, 
  rows_to_pack = rows_to_pack,
  file.path("Output", paste0("AllModels", suffix_onlypushing, ".xlsx")), 
  merge_option = 'both', 
  simplify_2nd_row = TRUE,
  line_above_rows = c(1,2),
  line_below_rows = c(-1)
)

  
summary_all_models
Subjective MVPA Hurdle Lognormal
Device-Based MVPA Log (Gaussian)
Mood Gaussian
Reactance Dichotome
Hurdle Component
Non-Zero Component
exp(Est.)_hu pa_sub 95% CI_hu pa_sub exp(Est.)_nonzero pa_sub 95% CI_nonzero pa_sub exp(Est.) pa_obj_log 95% CI pa_obj_log Est. mood_gauss 95% CI mood_gauss OR is_reactance 95% CI is_reactance
Intercept 3.79*** [ 2.28, 6.37] 57.44*** [49.33, 67.17] 124.09*** [109.75, 140.55] 3.77*** [ 3.54, 4.00] 0.40*** [0.22, 0.67]
Within-Person Effects
Daily individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA
Daily partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA
Daily individual’s experienced pushing 1.55** [ 1.17, 2.18] 0.97 [ 0.91, 1.03] 1.02 [ 0.97, 1.08] 0.02 [-0.05, 0.09] 1.36* [1.06, 1.84]
Daily partner’s experienced pushing 1.84*** [ 1.36, 2.72] 0.97 [ 0.90, 1.03] 1.02 [ 0.97, 1.06] 0.09* [ 0.00, 0.17] 1.15 [0.79, 2.07]
Day 0.61* [ 0.37, 0.99] 0.91 [ 0.79, 1.05] 0.97 [ 0.87, 1.08] 0.28** [ 0.11, 0.45] 1.78 [0.68, 4.81]
Own Actionplan NA NA NA NA NA NA NA NA NA NA
Partner Actionplan NA NA NA NA NA NA NA NA NA NA
Daily weartime NA NA NA NA 1.00** [ 1.00, 1.00] NA NA NA NA
Between-Person Effects
Mean individual’s experienced persuasion NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced persuasion NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pressure NA NA NA NA NA NA NA NA NA NA
Mean partner’s experienced pressure NA NA NA NA NA NA NA NA NA NA
Mean individual’s experienced pushing 0.47** [ 0.25, 0.80] 1.32* [ 1.01, 1.72] 0.97 [ 0.85, 1.10] 0.13 [-0.08, 0.33] 5.68*** [2.38, 16.59]
Mean partner’s experienced pushing 0.63 [ 0.34, 1.08] 1.17 [ 0.89, 1.53] 1.04 [ 0.92, 1.18] 0.09 [-0.11, 0.29] 2.41 [0.94, 7.27]
Mean weartime NA NA NA NA 1.00 [ 1.00, 1.00] NA NA NA NA
Random Effects
sd(Intercept) 1.26 [0.92, 1.75] 0.37 [0.28, 0.50] 0.31 [0.24, 0.42] 0.61 [0.48, 0.80] 1.53 [1.04, 2.25]
sd(Daily individual’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced persuasion) NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pressure) NA NA NA NA NA NA NA NA NA NA
sd(Daily partner’s experienced pressure) NA NA NA NA NA NA NA NA NA NA
sd(Daily individual’s experienced pushing) 0.47 [0.07, 1.02] 0.06 [0.00, 0.16] 0.09 [0.02, 0.17] 0.07 [0.01, 0.16] 0.17 [0.01, 0.55]
sd(Daily partner’s experienced pushing) 0.51 [0.09, 1.10] 0.09 [0.01, 0.18] 0.03 [0.00, 0.10] 0.14 [0.06, 0.24] 0.41 [0.03, 1.12]
Additional Parameters
sigma NA NA 0.67 [0.64, 0.70] 0.54 [0.52, 0.57] 0.91 [0.87, 0.94] NA NA 
report::report_system()

Analyses were conducted using the R Statistical language (version 4.4.2; R Core Team, 2024) on Windows 11 x64 (build 26100)

report::report_packages()
  • beepr (version 2.0; Bååth R, 2024)
  • R.methodsS3 (version 1.8.2; Bengtsson H, 2003)
  • R.oo (version 1.27.0; Bengtsson H, 2003)
  • R.utils (version 2.12.3; Bengtsson H, 2023)
  • brms (version 2.22.0; Bürkner P, 2017)
  • digest (version 0.6.37; Eddelbuettel D, 2024)
  • Rcpp (version 1.0.13.1; Eddelbuettel D et al., 2024)
  • bayesplot (version 1.11.1; Gabry J, Mahr T, 2024)
  • lubridate (version 1.9.3; Grolemund G, Wickham H, 2011)
  • DHARMa (version 0.4.7; Hartig F, 2024)
  • wbCorr (version 0.1.22; Küng P, 2023)
  • see (version 0.9.0; Lüdecke D et al., 2021)
  • tibble (version 3.2.1; Müller K, Wickham H, 2023)
  • R (version 4.4.2; R Core Team, 2024)
  • openxlsx (version 4.2.7.1; Schauberger P, Walker A, 2024)
  • ggplot2 (version 3.5.1; Wickham H, 2016)
  • forcats (version 1.0.0; Wickham H, 2023)
  • stringr (version 1.5.1; Wickham H, 2023)
  • rvest (version 1.0.4; Wickham H, 2024)
  • tidyverse (version 2.0.0; Wickham H et al., 2019)
  • readxl (version 1.4.3; Wickham H, Bryan J, 2023)
  • dplyr (version 1.1.4; Wickham H et al., 2023)
  • purrr (version 1.0.2; Wickham H, Henry L, 2023)
  • readr (version 2.1.5; Wickham H et al., 2024)
  • xml2 (version 1.3.6; Wickham H et al., 2023)
  • tidyr (version 1.3.1; Wickham H et al., 2024)
  • knitr (version 1.49; Xie Y, 2024)
  • kableExtra (version 1.4.0; Zhu H, 2024)
report::cite_packages()